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SECTION  I 


INTRODUCTION 


In  recent  years  satellite  communications  systems  have  received 
considerable  attention  by  both  commercial  and  military  users.  Much 
of  this  attention  was  due  to  the  supposition  that  the  frequency  ranges 
used  were  insensitive  to  the  ionospheric  propagation  conditions  with  the 
possible  exception  of  absorption  from  nuclear  bursts.  This  postulated 
insensitivity  made  satellite  systems  very  attractive  compared  to  other 
systems  such  as  HF  which  rely  upon  the  ionosphere  for  their  operation. 

The  great  promise  of  satellite  systems,  however,  failed  to  materialize 
due  to  the  scintillation  problem.  Scintillations  are  large  fluctuations 
in  the  amplitude  and  phase  of  a received  signal  caused  by  fluctuations 
in  the  ionospheric  electron  density  through  which  the  signal  passes. 

An  added  dimension  has  been  added  to  the  scintillation  problem  by  the 
realization  that  similar  propagation  conditions  might  result  from 
atmospheric  nuclear  bursts. 

These  scintillation  problems  have  created  concern  on  behalf 
of  the  military  services  as  satellite  communications  systems  are 
expected  to  take  an  ever  increasing  roll  in  military  communications. 

This  report  is  a result  of  that  concern.  Contained  herein  is  an  analysis 
of  two  communication  links  in  two  disturbed  ambient  environments,  along 
with  some  general  discussion  of  ionospheric  environments. 


The  two  systems  treated  are  an  M-ary  frequency  shift  key  system  (MFSK) 


t • 


and  a binary  phase  shift  system  (BFSK).  It  is  hoped  that  the  discussion 
and  calculation  results  will  be  useful  to  others  working  in  this  area  as 
well  as  better  defining  the  scintillation  difficulties  encountered  by 
these  systems. 

This  document  consists  of  seven  sections  excluding  the  introduction. 
The  second  section  discusses  the  propagation  methods  used  and  their 
justification.  The  third  and  fourth  sections  discuss  the  frequency  shift 
key  and  the  binary  phase  shift  key  models,  respectively.  The  determination 
and  a discourse  on  the  environments  used  are  found  in  the  fifth  section. 
Next  the  results  of  some  calculations  are  presented.  A brief  summary 
concludes  this  paper. 
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SECTION  II 


PROPAGATION 


With  the  advent  of  satellite  communication  links  and  the  discovery 
of  scintillated  signal  conditions  under  the  equatorial  and  polar  ionospheres, 
it  has  become  necessary  to  extend  VHF,  UHF,  and  higher  frequency  propagation 
techniques  (ref  1 to  8).  These  methods  were  generally  based  on  a thin  screen 
approximation  geometrical  optics,  weak  scattering,  or  some  combination  of 
these  approximations.  The  method  presented  here  is  not  restricted  by  the 
above  approximations  and  also,  for  many  situations,  supplies  more  information 
on  the  statistical  properties  of  the  propagated  wave  than  do  the  previous 


methods , 


There  is,  already  in  existence,  a method  by  Uscinski  (ref  9)  is 
not  subject  to  the  above  limitations.  This  previous  method  will  calculate 
the  same  information  as  the  method  presented  here,  but  appears  more 
difficult  numerically  to  evaluate. 

Let  the  signal  incident  on  the  ionosphere  be  represented  as 


E = (E  + E_  + i E )e 
R I 


where 


E = mean  field 

E^  = in  phase  quadrature  component 
Ej.  = out  of  phase  quadrature  component 
K = 2tt/A 


Let  me  define  two  statistical  quantities  called  the  complex  autocorrelation 


function  and  the  asymmetry  autocorrelation  function,  respectively,  where 


an  overline, 


, represents  the  ensemble  average: 


G(X,p)  = E*(X,Z)  E(X,Z+p) 


K + ED(X,Z)E  (X,Z+p)  + Et(X,Z)Et(X,Z+p) 


F(x,p)  = E(X,Z)E (X,Z+p) 


« E + ER(X,Z)ER(X,Z+p)  - EI(X,Z)EI(X,Z+p) 


+ 2i  ER(X,Z)EI(X,Z+p) 


since  ER(X,Z)  = E-^X.Z)  = 0 and  ER(X,Z)E];(X,Z+p)  = Ej.  (X,Z)ER(X,Z+p)  . 

The  symmetry  in  the  cross  correlation  of  ER  and  Ej.  will  result  from 
certain  assumptions  about  the  striated  medium.  The  second  transverse 
coordinate  will  be  suppressed  for  convenience.  Its  Inclusion  in  the 
expressions  that  follov?  is  a minor  matter. 

Tatarskii  (ref  10)  explicitly  derived  th  equations  for  the  mean  field 
and  the  complex  autocorrelation  function.  The  equation  for  the  asymmetry 
autocorrelation  function  can  be  derived  similarly  to  the  complex 
autocorrelation  function.  The  equations  are 


— 1 = - K2E(X)A(X,0)/2 


dF(X.p) 

dX 


A *ll£i£l  _ k2F(X,p)[A(X,0)  + A(X,p)  ] 
dp 


dG (X,p) 


- G(X,p)[A(X,0)  - A(X, p) ] 


where 


A(X,p)  = n(X,Z)n(X\Z+p)  dX' 

*'—00 

n(X,Z)  ■ index  of  refraction  fluctuations 

"n(X,Z)"  is  assumed  to  be  a zero  mean  gaussian  random  variable. 

These  conditions  on  n(X,Z)  are  sufficient  but  not  in  all  cases  necessary 

to  the  derivation  of  equations  5 and  6.  For  example,  if  a signal  suffers 

many  scatters  traversing  a striated  medium,  then  equations  5 and  6 can 

be  justified  by  using  the  central  limit  theorem.  In  any  event  experimental 
2 

work  at  SRI  has  shown  the  assumption  of  Gaussian  distributed  fluctuations 

2 

to  be  justified.  Equation  4 is  not  necessary  as  G(X,°°)  # E(X)  , but  it 
is  carried  in  calculations  to  examine  the  degree  to  which  the  boundary 
conditions  on  F(X,p)  and  G(X,p)  are  satisfied. 

The  limitations  on  the  solutions  of  equations  5,  6,  and  7 are 
discussed  in  Tatarskii  (ref  10)  and  will  be  summarized  here. 


X«i 

(8) 

o 

L>>£ 

(9) 

o 

X a<<l 

(10) 

°e2<<1 

(ID 

aLaQ2«l 

(12) 

L«X 

m 

(13) 
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where 


£ = striatlon  scale  size 

o 

L **  striatlon  layer  thickness 


a = mean  field  attenuation  coefficient 
0^2  = mean  square  scattering  angle 
X^  = characteristic  distance  for  attenuation 
by  back  scattering 

All  of  these  conditions  are  easily  met  for  present  problems  of 
interest.  An  implicit  ^sumption  in  using  the  A(X,p)  function  with  an 
explicit  X dependence,  is  that  the  variation  of  A along  X is  small  over 
a scale  size.  If  die  striated  medium  is  homogeneous,  then  A(X,p)  is 
symmetric  with  respect  to  P.  The  symmetry  in  P results  in  like  symmetry 
in  F(X, p)  and  G(X,P). 

The  functions  F(X,p)  and  G(X,p)  are  not  the  most  useful  qualities. 
Instead,  the  correlation  functions  of  the  field  components  are  the  ultimate 
product. 


ER(X,Z)ER(X,Z+p)  = [G(X,  P)  + Real(F(X,p))]/2  - E" 


EI(X,Z)EI(X,Z+p)  = [G(X,  p)  - Real(F(.X,p))]/2 


ER(X,Z)EI(X,Z+p)  = Imag(F(X,p))/2 


The  correlation  function  of  the  index  of  refraction  fluctuations 


used  in  all  present  calculations  is  of  the  form 


n(X,Y,Z)n  (X  + e,  Y + p,  Z + T})  = 


C (u)  n„(X)  6n  exp 


•# » ,r  wi 


where 


C(w) 


\mW  / 


np(X)  = mean  electron  density 


5,  = outer  scale  size 

o 


2 — 2 _ 2 
6 = local  value  of  (n  - n„  )/n 

n e e ' ' e 


The  one  sided  power  spectral  density  (PSD)  of  the  electron  density 
fluctuations  corresponding  to  equation  17  is 


PSD(K) 


2 2 3 

= ne(X)2  6/  *0 


2 2 2 2 
2tt2(1  + ) 


(1 


In  actual  ionospheric  experiments  the  electron  density  fluctuations  are 
measured  along  a given  path.  The  one -dimensional  form  of  (18)  is 


PSD(k^)  = 5n  ^o 


(1' 


t r(l/£0  + Kx  ) 


where  is  perpendicular  to  the  local  magnetic  field  and  K parallel  is 


assumed  to  be  zero.  This  last  assumption  is  often  very  reasonable  as 
the  ionospheric  structure  tends  to  elongate  along  the  ambient  magnetic 
field.  This  form  of  the  PSD  closely  fits  data  measured  by  several  authors 
(ref  11  *) . The  correlation  function  assuming  the  propagation  path  is 
perpendicular  to  the  magnetic  field  is  reduced  to 


•2  „ 2 


n(X,Y)  n(X  + e,  Y + P)  = C(w)  ne(X)  6n 


exp  [-(e2  + p2)1^2/£  ] . 


(20 


*Private  communication  with  M.C.  Kelley  and  K.C.  Yeh. 


This  form  is  used  in  equation  7 for  propagation  cclculations  near 
the  equator.  When  propagation  in  some  polar  scenarios  is  calculated,  the 
lines  of  sight  are  sometimes  more  parallel  than  perpendicular  to  the  field 


lines.  In  this  case  equation  17  is  used  with  x2  " P2  + n2,  and  the  transverse 
operator  in  equation  5 is  replaced  as  follows 

dp2  dx2  XdX  (21) 

In  the  analysis  of  the  communication  links  to  follow,  it  is  assumed 
that  Er  and  Ej.  are  Gaussian  random  variables,  TMs  assumption  is  reasonable 
when  the  propagation  disturbance  is  small  sit.  v<?  Ihe  propagation  is  still 
linear  and  Gaussian  distributed  electron  density  fluctuations  were 
assumed.  The  Gaussian  assumption  also  is  valid  for  the  strong  multiscatter 
case  as  the  resultant  propagated  signal  is  the  sum  of  many  independent 


scattered  signals.  In  between  these  two  extremes,  the  matter  is  not 
so  clear. 

Fortunately  the  fourth  order  moments  of  the  propagated  signals  for 
some  ionospheric  cases  have  been  calculated  by  Yeh*.  The  root  mean 
square  of  the  amplitude  fluctuations  was  deemed  to  be  the  most  useful 
measure  of  the  Gaussian  assumption.  The  maximum  error  in  the  mean 
square  amplitude  fluctuation  between  100  and  1000  MHz  was  about  18 
percent  with  the  value  using  the  Gaussian  assumption  being  the  largest. 
Fremouw  and  Rhino  (ref  12)  have  also  shown  that  scintillation  data  are 
consistent  with  the  Gaussian  assumption. 

It  is  important  to  remember  at  this  point  that  the  propagation  does 
calculate  the  first  and  second  moments  of  the  field  correctly,  and  the 


gaussian  assumption  on  E and  E is  an  added  feature  independent  of  the 

R 1 

propagation. 

All  of  the  propagation  discussion  thus  far  has  dealt  with  the 
propagation  of  a single  frequency.  With  the  consideration  of  spread 
spectrum  techniques,  it  is  necessary  to  know  the  correlation  of  signals 


of  different  frequencies  to  evaluate  these  techniques.  The  equation 
to  calculate  two  frequency  correlations  is  a generalization  of  equation 
6 and  was  first  used  extensively  by  Liu  (ref  13)  and  Ulaszek  (ref  14). 


H(X,p) 

dX 


= l(Kl~K2)  d2H(X,p) 


2KjK2 


l / («iW)  A (X.O)  _ A (X,p)  j H(X, 

(1K2  \ 2K1K2  / 


P)  = 0 


where 


H(X,p)  = E„  (X,Z)  E„  (X,Z+p) 


An(X,p)  = J 


Ane(X,Z)Ane(X+n,  Z+p)  dn 


Ane  (X,Z)  = electron  density  fluctuations 


, 2 _ 47ie2 
^ me2 


The  real  part  of  H(X,0)  is  taken  as  the  correlation  coefficient 
where  initially 


E (0,Z)  E (0,  Z+p)  = 1. 
K1  K2 


The  imaginary  part  is  the  sine  of  the  phase  difference  between  the 


two  frequencies  not  including  dispersion. 

In  real  practical  problems  the  electromagnetic  waves  are  not  plane 


waves  but  have  some  curvature.  Because  of  this,  it  is  necessary  to 


examine  the  limits  of  the  plane  wave  formalism.  The  following  will 
be  done  in  two  dimensions  to  illustrate  the  desired  limits.  The  general 
form  of  a cylindrical  wave  is 


E 


yQ(x,p) 


i K(X2+p2)li/2 
e 


(x2+p2> 


1/4 


(24) 


where  in  the  absence  of  propagation  effects 


h0(X,p)  = 1 


Equation  14  can  be  expanded 


E = Un(X,p)  e_ 


i 

2X 


Tfl 


iKX  iKX 

A = y(X,p)e 


X 


(25) 


(26) 


For  p _<  7 km  and  X = 250  km,  the  amplitude  of  the  approximation  is  in 
error  about  5 percent.  For  any  frequency  and  the  same  p and  X,  the 
radius  of  the  wave  front  is  off  by  £ 4 cm.  Thi9  error  is  insignificant 
for  present  purposes,  even  though  this  may  be  several  radians  in  phase 
at  higher  frequencies.  This  is  because  we  are  interested  in  the  incoherent 
scattering  of  energy  where  kilometers  of  propagation  path  are  needed  to 

see  any  effects.  The  complex  correlation  function  is 

i_K  (p  2 - p 2) 

2X 


* 

y 


(x,p2)y(x,p1) 


yQ(x,p2)ty:x,p1) 


x 


If  z = P1  - 


P2  and  R =4y  (Px  + P2>, 


y (X,p2)y(X,p1)  = G(X,R,Z) , and 


(27) 


i KRZ 

X 


Mo(X,P2)yo(X,Pl) 


G (X,R,Z) , then  G(X,R,Z) 
o 


G (X,R,Z)- 
o 


(28) 


From 


Tatarskii  (ref  10),  the  equation  for  G(X,R,Z)  is 


^^-|IzIrg(x‘r’z) 


+ K2  [A(X,0)  - A(X,Z) ] G(X,R,Z)  = 0 
imit  of  a plane  wave  3G(X,R,Z)  =0,  and  equation  6 is  obtained. 


(29) 


In  the  limit 

OR 

For  initial  conditions  on  C(X,R,Z) , C0<X,R,Z>  is  sot  to  one.  It  is 
easy  to  show  that  if  there  is  no  scattering,  that  is,  A(X,Z)  - 0.  then 
equation  28  is  a solution  of  equation  29.  This  results  because  equation  26 
is  an  exact  solution  to  the  parabolic  equation. 

Substituting  equation  28  into  equation  29, 


dGo(X,R,Z)  _ i 

dx K 


a2G0(X,R,Z) 

3Z3R 


iK  /r  3 Gn(X,R,Z)  ^ Z 9Gp(X,R,Z)  j 
iT  \ 3R  3Z  L 


+ K2[A(X,0)  - A(X,Z) ] G0(X,R,Z)  = 0 

Since  the  initial  condition  on  Gq  has  no  R dependence,  it  can  never 
develop  any . Thus 


(30) 


dG0(X,R,Z)  + z 3Gq(X,R,Z) 
dX  ’ x 3Z 


+ K2  [A(X,0)  - A(X,Z) ] Gq(X,R,Z)  = 0 (31) 


11 





G (X,R,Z)  at  large  Z is  the  fraction  of  power  not  scattered  and  is 


identical  with  the  plane  wave  result.  The  second  term  in  equation  31 
is  a geometrical  divergence  term  . Thus  the  solution  to  equation  31 


looks  approximately  like  the  plane  wave  solution  with  the  transverse 


Xo 


coordinate  scaled  as  -tt—  where  Xn  is  the  center  of  the  scattering 


layer. 


The  asymmetry  auto -correlation  function  is  . .2  2. 

i +p2  ) 


2X 


y(X,p2)y(X,pi)  = h0(X,p2)yo(X,p1) 


(32) 


Thus,  as  before, 


i K (2R2+Z2/2) 


2X 


F(X, R,Z)  = F (X,R,Z)  ±_ 


(33) 


The  equation  for  F(X,R,Z)  is 


3F(X,R,Z)  1 

3X  2iK 


2 32  , 1 32 

L 3Z2  2 3R2. 


F(X,R,Z) 


+ YT  [ A(X,0)  + A(X,Z) ]F(X,R,Z)  = 0 


(34) 


As  before  if  A(X,Z)  = 0,  then  F(X,R,Z)  is  an  exact  solution  to  equation  34 
with  FQ(X,R,Z)  = 1 as  the  boundary  condition.  The  equation  for  F (X,R,Z) 


is 


3F0(X,R,Z)  _ i 32Fo(X,R,Z)  z 3Fd(X,R,Z) 

ktz rr  v • 


8X 


K 3Z2 


3Z 


+ K [A(X,0)  + A(X,Z)]Fd(X,R,Z)  = 0 


(35) 
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Fq(X,R,Z)  is  not  a function  of  R since  A(X,Z)  is  nof.  Again,  there  is  the 

plane  wave  equation  plus  a geometric  divergence  term.  As  in  the  complex 

autocorrelation  function,  Fo(X,R,Z)  looks  very  much  like  the  plane  wave 

solution  with  the  transverse  coordinate  scaled  as  Xq/X.  Thus  the  pictures 

for  F (X,R,Z)  and  G (X,R,Z)  are  the  same, 
o o 

For  spherical  waves  with  an  axisymmetric  electron  fluctuation 
correlation  function,  equations  31  and  35  remain  the  same  except  that 

g2  ^ g 

— j goes  to + — irr  and  that  R and  Z are  now  radial  coordinates. 

3Z  3Z2  L dL 

To  visualize  the  final  physical  result,  imagine  the  solutions  to  equations 
31  and  35  being  wrapped  along  the  curved  wave  surface.  This  wrapping 
is  accomplished  mathematically  by  the  phase  terms  in  equations  28  and  33 
The  conclusion  of  this  analysis  is  that  the  plane  wave  results  with 
Xq/X  scaling  varies  little  from  the  cylindrical  and  spherical  wave  results. 
Calculations  verify  this  result. 

In  summary  adequate  methods  exist  to  calculate  the  propagation  effects 
of  a turbulent  ionosphere  on  traversing  signals.  The  included  equations 
allow  the  calculation  of  the  first  and  second  order  signal  statistics  at 
one  frequency  as  well  as  two  frequency  statistics.  Fourth  order  moments 
of  the  signal  are  calculable  as  in  the  calculations  of  Yeh*.  In  general, 
however,  the  calculation  of  the  fourth  moments  are  very  expensive;  thus 
the  use  of  the  Gaussian  assumption  on  the  quadrature  components  of  the  signals. 

The  calculation  methods  described  here  are  relatively  inexpensive. 

In  general  less  than  3 minutes  of  CPU  time  on  a CDC  6600  were  needed  per 
calculation  for  either  plane,  cylindrical,  or  spherical  geometries. 

*Private  communication  with  K.C.  Yeh. 
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M-ARY  FREQUENCY  SHIFT  KEY  DEMODULATOR 


The  received  signal  is  assumed  to  be  of  the  form 

i2ir(f  c + + fQ)t 

S (t)  = [E  + ER(t)  + iEj-(t)  ]e 
n 


+ N(t) 


C36) 


where  E = mean  field 

E (t)  = in-phase  component 
R 

= out— of— phase  component 

f = carrier  frequency 
c 

Af  = frequency  shift 
f = carrier  frequency  error 
n = frequency  shift  number 


N(t)  = white  gaussian  noise 


F„r  a binary  system  n equals  ± 1,  and  for  an  8-ary  system  n equals 
+1,  +3,  +5  and  ±7.  The  in-phase  and  out-of-phase  components  are  assumed 

random  Gaussian  variables. 

The  demodulator  is  a group  of  matched  filter-rectifier  pairs  in 

parallel.  Figure  1 contains  a binary  system.  The  operation  of  the 

matched  filter  matched  to  the  incoming  signal  is 

— . / r _i_  m J 


Rn 


= Eo  f 

N J 


-i  2ir(fu  + n_M)t 

S(t)  i dt 


(37) 


= dumped  filter  output  at  T 
N = unilateral  noise  power  spectral  density 
E = expected  signal  strength 


E0  rT  -i  2n(f  + M.)t  Ri 

~ / S(t)e  ' C 2'  dt  


Ec  fT  -1  2*1 

f j S(t,e 
0 

<t  -Mk 

“ 2'dt 

R1  . 

IR1  1 

Ip  1 

1 r2  1 

lRi I $ M 


Figure  1 Binary  FSK  demodular 

The  operation  of  the  rectifier  following  the  matched  filter  is 
represented  by  taking  the  absolute  value  of  R^.  The  output  of  the  m 
filter,  assuming  the  n signal  is  being  transmitted,  is 


E„rTr.  ......  ,.i  i!,M*  £>. 


E + FR(t)  + iEj(t) J e 


E ' rT  -i  27i/f  +m^!)t 

+ r~  I N(t j e ' c 2/  dt  (38) 

•’o 

The  outputs  of  the  filters  consist  of  the  sum  of  two  complex  numbers. 
The  first  is  a deterministic  term  dependent  on  the  received  signal,  the 
particular  filter,  and  the  total  phase  error.  The  second  term  is  the 
random  noise  term.  Figure  2 portrays  an  example  of  the  output  of  the 
filters  of  a binary  system  with  n equal  1.  The  vectors  from  the  origin, 
r^  and  ^ represent  the  deterministic  terms,  and  the  second  vectors  are 
the  noise  terms. 

From  the  previous  discussion  it  is  clear  that  the  noise  vector  is 
random  in  orientation  with  a Rayleigh  distributed  amplitude.  Given 
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Figure  2.  Sample  Filter  Output 


r , r2>  and  knowing  the  noise  statistics  the  probability  of  error,  i.e. 
| R2 | > |Rx| , can  be  calculated  (ref  15).  Let 

d2  = E 2T 
o 

N 


fl  d 

e ' l~2l 

fl  d 

Then  the  probability  of  error  is 

2 2 

Pe  = 1 - Q(a,3)  +~e~  &+LL  Iq  (op)  1 

~(q2+g2)  f g /g\2  1 

1-Q(a,g)  - e 2 j^a  Ij(oi3)  ^2  + . . .J  < 

"d2"  is  the  signal  to  noise  ratio.  It  is  also  twice  the  baud  energy-to- 
noise  ratio..  In  (X)  is  the  conventional  integer  order  modified  Bessel 


function. 


The  simulation  of  an  M-ary  system  proceeds  as  follows.  First 
using  the  method  described  in  appendix  A,  a random  sample  of  E + E^(t) 

+ iEj(t)  is  generated.  Given  f and  n designating  the  signal  present, 

R is  determined  at  appropriate  sampling  times.  For  each  sampling  time  the  prob 
ability  of  error  is  one  minus  the  product  over  m ^ n of  the  probabilities  that 

J Kni  | < | R [ . This  probability  of  error  at  each  sampling  time  is  then 
used  to  determine  various  performance  parameters  of  the  system.  The 
channel  containing  the  signal  is  chosen  randomly.  Meaningful  statistics 
are  achieved  by  processing  a sufficient  number  of  sample  signals. 

This  methodology  allows  simulation  of  all  the  major  aspects  of  this 
type  of  receiver.  In  particular  nonorthogonality  of  channels  and 
frequence  error  in  the  incoming  carrier  are  included.  The  simulations 
are  inexpensive,  consuming  less  than  3 minutes  of  CPU  times  on  a CDC  6600. 
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SECTION  IV 


BINARY  PHASE  SHIFT  KEY  DEMODULATOR 


The  received  signal  is  represented  as 


E ' (t)  = M(t)  [E  + ER(t)  + i Ej(t) ) e 


i 27Tfct 


(44) 


where 


E = mean  field 


E (t)  = in  phase  component 
R 


Ej(t)  = out  of  phase  component 


f = carrier  frequency 
c 


M(t)  = modulation 


M(t)  is  either  +1  during  each  bit  period  corresponding  to  a 0 or  a 
1 bit.  For  purposes  of  analysis  only,  the  real  part  of  equation  44  is 
necessary.  E (t)  and  E-^t)  are  assumed  Gaussian  distributed.  Let 


s(t)  = Real  [E  (t) ] 

= M(t)  E(t)  cos  (8(t)  + 27Tfct) 


(45) 

(46) 


where 


E(t)  = [<E  + ER(t))2  + Ej2 (t) ] 1/2 


0 ( t ) = tan"1[Ei(t)/(E  + ER(t) ) ] 


The  model  demodulator  is  described  in  figure  3. 
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Figure  3.  Demodulator  Model 

The  incoming  signal  is  first  filtered  by  an  IF  filter  with  a one  sided 
bandwidth  of  W.  The  signal  then  goes  to  a multiplier  and  a doubler  phase- 
lock  ed-loop  (DPLL) . The  DPLL  generates  a coherent  reference  which  is  fed 
into  the  multiplier  and  is  multiplied  by  the  original  signal.  The 
multiplier  result  is  integrated  over  a baud  period  and  the  data  determined. 
This  is  essentially  the  same  model  as  in  reference  9. 

The  demodulator  description  will  consist  of  two  sections.  The  first 
will  discuss  the  IF  filter  and  the  doubler  phase-locked-loop . The  second 
will  treat  the  multiplier  and  the  integrator.  The  initial  discussion  to 
follow  parallels  closely  the  discussion  in  reference  9.- 

IF  Filter  and  DPLL 

The  IF  filter  is  assumed  to  have  the  following  transfer  function. 
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where  fQ  is  the  carrier  frequency  and  jp-J  i-8  shown  in  figure  4« 


Figure  4^  tt  Function 

The  input  noise  spectrum  is  assumed  constant  over  the  bandwidth  of 
the  filter  thus  the  power  spectral  density  of  the  filtered  noise  is 


where  Nq  is  the  unilateral  spectral  density  of  the  noise 

The  first  element  in  the  DPLL  is  a square  law  device.  This  device 
removes  the  BPSK  modulation.  With  an  input  signal  of 

v(t)  = s ( t ) + n(t) 


s(t)  = M(t)  E(t)  cos  (0 (t)  + 2ir  f t) 


the  output  is 


2,  x E2(t) 

s^(t)  =2  [i  + cos  (4tt  fQt  + 26(t))  ] 

+ 2M(t)  E(t)  cos  (27TfQt  + 9(t))n(t)  + n2(t) 

The  last  two  terms  in  equation  51  are  noise  derived  terms, 
power  spectral  density  of  these  terms  is  derived  in  appendix  B. 

Sn2(f)  = 4Sn(f)  * Sg(f)  + 2Sn(f)  * Sn(f) 

+ DC  terms 

where  Sg(f)  is  the  power  spectral  density  of  the  signal,  s(t). 
sn(f)  * Sn(f)  term  is  readily  evaluated. 

2 

2Sn(f)  * Sn(f>  - !!o_«  L ^£-2£o  j 

+ M^)  + 2A(ff)] 

where  the  A function  is  illustrated  in  figure  5, 


tAmJk 


Let  us  now  look  at  the  S x S term  where 

• n s 


s(t)  = M(t)  E(t)  cos  (2irf0t  + 0(t)) 

If  E(t)  and  0 ( t ) are  assumed  to  change  slowly  with  respect  to  the 
carrier  frequency,  then  S„(f)  can  be  written  as 

O 

Vf>  ■ i tV£-£o>  + sb'-£-£o>) 

where  Sb(f)  = TT  ^ Sc(f) 


and  S (f)  is  the  power  spectral  density  of  M(t)E(t)e^^ . As  W decreases 
c 

some  of  the  signal  power  is  lost.  Let  the  fraction  of  signal  power  that 
is  lost  be  . Then 


sb(f)df 


N£/ 


S0(f)df 


where  \ Sc(f)df  is  twice  the  received  signal  power. 


From  equations  48  and  55  the  S x S term  can  be  evaluated. 

ns 


where 


Ss(f)  * Sn(f)  = g-  [Sbn(f+2f0) 


+ Sbn(f-2f0)  + 2Sbn(f)] 


sbn(f)  =sb(f)  **(§) 

W/2 

= j Sb(f-y)dy 


-W/2 


1 


From  equations  53  and  58  the  total  power  spectral  density  of  the 


noise  out  of  the  square  law  device  is 

N 

S„2(f>  ' ~ [Sbn(£+2fi)  + Sbn(£-2£o)] 


+ low  frequency  terms 


(60) 


The  total  noise  power  out  of  the  square  law  device  in  the  side 
bands  is 


/Sn2(f)df  = 2N02W2  + 2N0WQ(t)Nf 


(61) 


where  Q(t)  is  the  received  signal  power.  The  time  dependence  in  Q(t) 
denotes  the  variation  of  signal  power  from  scintillations,  absorption, 
and  such  effects.  In  terms  of  E(t) 


Q(t)«E2(t)/2  (62) 

In  the  discussion  of  the  phase  lock  loop  (PLL)  to  follow,  it  will  be 
noted  that  only  the  input  frequencies  in  a narrow  band  B^  centered  about 
+ 2fQ  are  important.  is  called  the  loop  bandwidth.  Since  is  usually 
much  smaller  than  W,  ■ the  noise  output  out  of  the  square  law  device  is 
characterized  as  gaussian  with  a unilateral  power  spectral  density  of 

N = 2N0Q(t)  Ng  + N02W  (63) 
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The  x noise  output  is  not  actually  gaussian  but  has  a x 
distribution.  If  < W,  however,  the  terra  can  be  considered  gaussian 
since  the  PLL  operation  is  equivalent  to  an  integration  over  many 

decorrelation  times  of  the  noise  Input.  The  output  signal  of  the  square 
law  device  is  now  represented  as 


S2(t)  = /2P(t)  sin  (2irf  t + <Kt)) 

C 0 

+ n’(t) 

where  P(t)  = E^*(t)/8 

f = 2f 
c o 

4>(t)  = 20 (t)  - tt/2 

n'(t)  = noise  characterized  by  equations  60  and  63 
The  change  in  symbols  is  to  conform. to  the  discussion  in  reference  20 
where  the  phase  lock  loops  of  present  concern  are  described  in  greater 
detail.  Figure  6 is  a diagram  of  the  phase  lock  loop  to  be  modeled. 


/P(t)  sin  (4) (t)  -$(t))  + n" (t)  + high  frequency  terms 


Figure  6 . Phase  Lock  Loop 


n"(t)  is  a zero  mean  Gaussian  process  with  bandwidth  W,  spectral  height 

A 

N/2,  and  center  frequency  zero.  In  figure  6,  4>(t)  is  the  phase  estimate 
of  the  incoming  signal.  The  high  frequency  terms  out  of  the  multiplier 
will  be  filtered  out  so  they  will  be  dropped.  The  output  of  the  voltage 
controlled  oscillator  (VCO)  is  a sinusoid  with  a frequency  shift  proportional 
to  the  input . The  output  is  written  as 


•Jl  cos  (27Tf  ct  + a f(u) 


= cos  (2irf  t + <f>(t)) 


where  f(u)  is  the  input  to  the  VCO.  "a"  will  be  taken  as  equal  to  one. 

Two  filters  will  be  considered.  The  first  is  just  a gain,  K /✓f-,  where 
PQ  is  the  value  of  P(t)  with  no  propagation  effects.  The  analysis  of 
this  first  order  loop  is  straightforward 

f (t)  = K^PU7  sin  [<Kt)  " i(t)]  + n"(t) 


= d<j>(t) 

dt 


Let  e(t)  = (J)(t)  - <j>(t) 


dSiEi  - sin  [e(t>]  - *=r  „»(t)  + 


The  second  filter  has  a response  of  — ^1  + ~ 


This  is  known  as  a second  order  loop. 

The  last  operation  is  a frequency  divider  which  results  finally  in 
the  desired  reference  signal,  r(t). 

I * • r(t)  = /2  cos  ^ 27TfQt  + 0 ( t ) - e(^-)  (69) 

1 

The  error  in  the  reference  generated  by  the  DPLL  is  thus  found  by  solving 
for  e(t).  The  solution  techniques  used  to  solve  equations  24  and  25 

I 

are  described  in  appendix  C. 
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Multiplication  and  Integration 

In  the  last  section,  solving  for  the  coherent  reference  r(t)  was 
described.  Multiplying  the  incoming  signal  and  integrating  over  the  baud 


period,  T,  results  in  the  statistic  R. 

T. 

R = 2 Jll(t)  E(t)  cos  ( 2tv f 0 1 + 0(t))  cos  ( 2iTf Dt  + 6(t)  - ) dt 

0 


+ 2 J n(t)  cos  (27TfQt  + 6 ( t ) - dt 


(70) 


If  R > 0,  then  the  phase  of  E(t)  is  taken  as  0 . If  R < 9 then  the 
assumed  phase  is  it.  The  last  integral  on  the  right  in  equation  28 
is  a random  zero  mean  Gaussian  variable  with  variance  NT.  The  decision 
process  can  be  written  as'  follows. 


d + G(0)  = R/(NT)1/2  i 0 


(71) 


where  d = ^ | E(t)  | cosj^e(t)Jdt/(NT) 


1/2 


G(0)  is  a Gaussian  sample  of  unit  variance.  Pe>  the  bit  error  probability, 
is  now  easily  found. 


Pe  = erf c (d)  = 


/2F  J 


- x2/2 


dx 


(72) 
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The  general  method  of  analysis  goes  as  follows.  Using  methods 


described  in  appendix  A random  ensembles  of  the  input  signals  are 
generated  and  $(t)  is  found  for  the  sample.  Next  the  appropriate  loop 
equation  is  solved  to  get  the  phase  error  as  a function  of  time.  Lastly 
P£  is  determined  at  appropriate  sampling  times  and  statistics  on  the 
error  performance  of  the  system  are  performed.  Adequate  statistics 
are  insured  by  generating  a sufficiently  large  number  of  sample  signals. 

This  model  provides  an  inexpensive  simulation'  of  all  of  the  major 
aspects  or  the  given  BPSK  demodulator  including  the  effects  of  noise 

and  nonlinearities  in  the  PLL.  The  simulation^  required  less  than  8 

/ 

minutes  of  CDC  6600  CPU  time  for  each  case.  / 


SECTION  V 

PROPAGATION  ENVIRONMENTS 


i 


I 

1 


Global  morphology  studies  of  scintillations  at  AFCRL  (ref  18,19,20,21) 
indicate  that  scintillations  are  predominately  equatorial  and  polar 
phenomena.  The  environments  used  in  this  study  are  intended  to  be 
representative  of  these  two  regions.  The  equatorial  scintillations 
predominately  occur  within  15  degrees  of  the  magnetic  equator  between 
1900  and  0100  local  time.  They  are  present  year  around  with 
maximums  in  frequency  near  the  equinoxes.  Some  polar  scintillation 
is  virtually  always  present  above  70  degrees  geomagnetic  latitude 
which  is  the  scintillation  boundary  at  about  0900  local  time.  The 
boundary  moves  south  with  time  reaching  about  57  degrees  geomagnetic 
latitude  at  2100  hours.  The  boundary  is  generally  less  than  60  degrees 
between  1800  and  2400  hours  local.  During  disturbed  magnetic 
conditions  the  scintillation  boundaries  can  go  further  south.  The 
severity  of  scintillations  is  not  uniform  in  the  polar  region.  The 
most  severe  occurances  occur  during  the  late  evening  hours  in  a region 
bounded  on  the  south  by  the  already  mentioned  boundary  and  extending 
northward  5 to  10  degrees. 

In  section  II,  the  environmental  description  necessary  to  do 
the  propagation  section  included  the  following  parameters. 

ng(X)  = mean  electron  density 

(n  (X)^  - n (X)^n  (X)^  = a ^ = relative  fluctuation  power 
e e e n 

i = outer  scale  size 
o 


1 


The  mean  electron  density  and  the  relative  fluctuation  power  (RFP)  are 

treated  as  separate  parameters  primarily  because  the  available  in  situ 

2 . 
data  are  in  terms  of  0^  and  not  including  ng(X^ . 
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As  the  available  in  situ  data  was  examined  it  became  clear  that 


insufficient  data  existed  to  unambiguously  specify  the  propagation 

environments.  In  the  equatorial  region,  there  are  two  sets  of  data 

available.  The  first  set  by  Kelly  and  Mozer*  is  from  rockets  launched 

from  Natal,  Brazil,  and  contains  the  only  F region  peak  data  known  to  the 

2 

author.  From  these  data,  estimates  of  a /£  and  £ can  be  made.  The 

n o o 

2 

estimates  of  a /£  are  shown  in  table  1.  The  second  set  of  equatorial 
no 

data  is  from  Dyson  (ref  11)  and  is  from  446  km,  above  the  F region  peak. 

2 

Only  a value  for  a /£  can  be  estimated  and  it  is  also  in  table  1. 
n o 

More  satellite  and  rocket  in  situ  data  exist  but  thus  far  it  has  not 
yet  been  reduced. 

In  the  polar  regions  the  data  are  from  Sagalyn*  and  Dyson  (ref.  11). 

data  unfortunately  only  include  altitudes  above  about  700  km.  For  these 
2 

data  both  £ and  o /£  can  be  estimated.  As  seen  in  figure  7,  £ and 
o n o ° ’ o 

2 

/£q  are  larger  and  smaller  respectively  than  their  lower  altitude 
counterparts  near  the  equator. 

TABLE  I.  ENVIRONMENT  PAP \METERS  (EQUATORIAL) 

£ (Km)  a 2/£_(km 

o n 0 

Kelley  1.25  .126  < 0 2/£  < .66 

— no  — 

Dyson  > .45  .22 

Two  other  pieces  of  data  are  necessary  to  do  the  propagation  and 
the  subsequent  systems  analysis.  First  an  estimate  of  np(X),  the  mean 
electron  density  is  needed.  The  ionospheric  profiles  used  for  the 
equatorial  and  polar  cases  are  shown  in  figure  8.  It  is  also  necessary 
*Private  communication  with  M.C.  Kelly  and  R,  Sagalyn. 
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These 
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I 

I 

I 


i 


I 


I 1 


to  have  an  estimate  of  the  drift  velocity  of  the  ionospheric  structure. 
Measurements  by  Lincoln  Laboratories*  and  Paulson  and  Hopkins  (,ref  11)  near 
the  equator  indicate  that  150  + 75  meters/second  are  reasonable  estimates. 
These  values  are  also  typical  for  polar  regions*. 

2 

In  the  present  calculations,  due  to  the  paucity  of  data,  a /£  and  Z 

are  assumed  constant  for  the  entire  bottom  side  of  the  ionosphere. 

2 

This  use  of  the  constant  single  values  for  a /£  and  £ may  at  first 

° n o o J 

seem  a poor  representation  of  the  ionospheric  structure.  There  is  reason 

to  believe,  however,  that  single  values  are  not  unreasonable.  As  already 
2 

noted  £ and  0 /£  tend  to  be  larger  and  smaller,  respectively,  at  high 

o no 

latitudes  and  altitudes  than  the  corresponding  quantities  at  lower 

altitudes  and  latitudes.  From  general  considerations  of  the  gradient 

drift  and  other  instabilities,  the  author  believes  that  the  differences 
2 

in  £q  and  0n  /£q  are  primarily  altitude  related  since  both  these  quantities 
at  saturation  are  a function  of  the  ionospheric  gradients  giving  rise 
to  the  instability  (See  Sleeper,  ref  23).  The  steeper  gradients  appear  to 
peak  in  the  lower  F region  in  the  late  evenings  hours  in  both  the  polar  and 
equatorial  regions.  The  steep  gradients  are  most  likely  a result  of 
recombination  bottom  upward  of  the  F region  as  the  solar  ionization 
decreases.  (Note  that  no  mention  has  been  made  to  the  mechanisms  actually 
driving  the  instability.  This  is  because  the  saturation  of  the  instability 
is  not  dependent  on  the  driving  mechanism  but  only  the  initial  gradients 
and  the  mean  electron  density).  The  increased  RFP  and  the  maximum  electron 
density  of  the  F region  peak  conspire  to  make  this  region  dominant  in 
the  propagation  effects.  Thus  using  the  environment  parameters  representative 
of  the  F region  peak  is  reasonable.  This  postulated  dominance  of  the  F peak 
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in  producing  propagation  effects  has  been  questioned  by  some  (ref  2b)  but  all 

measurements  of  the  height  of  the  perturbing  regions  known  to  the  author 

indicate  an  altitude  of  about  300  km  or  less.  The  presence  of  Fresnel 

filtering  in  much  propagation  data  also  argues  against  thick  striated 

layers.  Hopefully,  existing  in  situ  data  for  equatorial  high  altitudes 

and  for  the  polar  F region  peak  in  the  F layer  irregularity  zone  will 

soon  be  reduced,  and  the  F layer  peak  question  can  be  resolved. 

2 

The  use  of  the  specified  values  of  an  /£q  for  the  polar  regions  is  of 
unknown  validity.  It  is  believed  by  the  author  that  the  numbers  are  reason- 
able for  the  polar  case  but  not  necessarily  optimum  because  the  detailed 
gradient  structure  in  the  polar  ionosphere  does  vary  significantly  from 
that  near  the  equator. 

2 

The  actual  determination  of  the  values  for  CT  /£  and  £ used  in 

n o o 

the  following  relevations  resulted  from  an  iteration  between  reasonable 
2 

values  for  a /£  , n (X) , and  £ and  between  extensive  propagation  data 
n o e o 

taken  by  Paulson  and  Hopkins  (ref  22,25,26)  at  selected  times  from  1970- 

1972  near  the  equator.  They  assumed  that  scintillation  existed  when  the 

fades  exceeded  6 dB  with  some  regularity.  Three  types  of  data  were 

used  here.  Cumulative  amplitude  distributions  for  250  MHz,  statistics 

on  fade  duration  for  250  MHz,  and  the  1 to  99  percent  points  of  some 

2 

2.3  gigahertz  peak  to  peak  data  were  used  to  estimate  a n /£q  and  £q 

2 

assuming  n (X)  of  figure  8.  It  should  be  noted  that  the  value  of  a 
e n 

can  be  varied  against  ng(X)  without  changing  the  net  propagation  effects. 

Thus  there  is  some  flexibility  in  picking  these  parameters.  In  addition 
the  choice  of  150  meters  per  second  for  the  ionospheric  drift  velocity 
influenced  the  final  numbers.  The  final  numbers  found  that  agree  with 
the  to  be  described  propagation  are 
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6 2/l  = 8.6%  km'1 
n o 

Z - 0.44  km 
o 

These  numbers  represent  a considerably  less  severe  propagation  environment 
than  that  represented  by  the  Kelley  data.  Considering  the  great  variability 
of  all  the  environment  parameters  in  both  time  and  space  and  the  fact 
that  the  propagation  data  are  averaged  through  and  over  different  nights, 
this  difference  is  probably  not  significant.  By  varying  the  ionospheric 
drift  velocity,  the  mean  ionospheric  profile,  and  the  regions  of  the 
ionosphere  that  are  unstable,  the  values  above  can  be  varied  over  considerable 
ranges  and  still  accurately  reproduce  the  propagation  data.  All  that 
can  be  said  is  that  no  gross  inconsistency  has  been  revealed  in  this 
analysis  between  the  propagation  and  environment  data.  Additional  data 
on  the  properties  of  the  ionosphere  must  be  accumulated  before  more  can 
be  said. 

Figure  9 contains  the  fade  duration  statistic  averaged  over  all  of 
the  Paulson  and  Hopkins  dtta  sets.  The  agreement  is  reasonable.  With  the 
given  data,  the  1-99  percent  peak-to-peak  amplitude  numbers  were  24  and 
2.6  dB  respectively  for  250  MHz  and  2.3  gigahertz,  respectively,  which,  is 
reasonable  agreement  with  Paulson  and  Hopkins.  Figure  10  is  a comparison 
of  experimental  and  calculated  cumulative  amplitude  data.  Finally  Paulson 
and  Hopkins  report  a correlation  of  0.6  for  two  frequencies  near  250  MHz 
separated  by  50  MHz.  The  calculated  value  was  0.7.  Thus  f-?r  the. 
equatorial  region  a consistent  and  reasonable  set  of  environmental 
descriptors  has  been  found.  The  use  of  data  from  two  widely  spaced 
frequencies  also  adds  confidence  in  using  those  numbers  to  calculate  for 
other  frequencies. 
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The  problem  of  finding  descriptors  for  the  polar  region  was  hindered 
by  a lack  of  data  either  in  situ  or  propagation.  As  already  mentioned 
the  in  situ  data  are  only  from  very  high  altitudes,  well  above  the  F region 
peak.  The  only  multifrequency  propagation  data  are  from  Pope  and  Fritz 
(ref  27).  These  data  was  taken  at  137  and  1695  MHz  by  examining  Automatic 
Gain  Control  (AGC)  records  for  ITOS  I and  ESSA  9 satellite  systems.  The 
137  MHz  data  are  Rayleigh  scattered  agreeing  with  my  calculations  using  the 
assumed  parameters.  The  1695  MHz  data  are  of  limited  use  here  for  several 
reasons.  The  most  important  is  that  the  S band  transmissions  were  made 
when  propagation  conditions  were  most  favorable  in  terms  of  lack  of 
scintillation.  Thus  these  data  are  probably  not  representative  of 
scintillated  conditions.  The  use  of  the  environmental  descriptors  from 
the  equator  in  the  polar  regions  is  not  inconsistent  with  these  data,  but 
that  is  all  that  can  be  said  until  additional  data  are  available. 


SECTION  VI 

SAMPLE  RESULTS 


To  illustrate  typical  calculations,  the  performance  of  a Binary 
Frequency  Shift  Key  (BFSK)  and  a Binary  Phase  Shift  Key  (BPSK)  system 
were  calculated  in  the  polar  and  equatorial  environments  previously 
described.  The  scenarios  are  a synchronous  satellite  down  to  an  aircraft 
going  300  m/sec  at  the  equator  and  in  the  polar  region  beneath  the 
north  wall  of  th  ionospheric  trough.  The  frequency  is  250  MHz.  Figure 
11  contains  the  calculated  correlation  functions  of  the  signals.  The 
equatorial  case  is  a Rayleigh  scattering  case.  The  polar  case  is 
neither  Rayleigh  or  Rician  since  the  variance  of  the  in  phase  and  out 
phase  components  are  not  equal. 

_2 

Both  modems  had  a 200  Hz  bandwidth  and  a bit  period  of  10 

seconds.  The  loop  bandwidth  of  the  PSK  receiver  is  12  hertz  (a=16, 

K=32)  . Figures  12  and  13  contain  the  results  of  the  link  calculations. 

-4 

Shown  are  the  signal  to  noise  ratios  necessary  for  a 10  bit  error 
rate,  the  bit  error  rates  for  all  the  cases,  the  maximum  error  free 
transmission  time  over  the  total  simulation  time,  and  the  window 
statistics.  The  meaning  of  the  window  statistics  cin  be  interpreted 
as  follows.  For  the  FSK,  polar,  and  15  dB  calculation,  the  point 
at  AT=5  seconds  means  that  in  a minute  there  will  be  three  intervals  of 
5 seconds  or  greater  with  no  bit  errors. 

At  15  dB  the  PSK  link  is  superior,  but  this  trend  is  reversed  at 
25  dB.  This  is  contrary  to  the  results  using  a linear  analysis  of  the 
PSK  receiver  illustrating  the  necessity  of  treating  the  nonlinearities 
of  the  PSK  receiver  with  noise.  A valuable  lesson  to  be  learned  from 


these  simple  results  is  that  even  though  the  bit  error  probabilities 
are  relatively  large,  significant  information  could  be  transmitted 
provided  the  message  is  appropriately  structured.  At  present  various 
schemes  are  being  considered.  One  possible  way  is  to  interleave  the 
bit  stream  over  a suitable  interval  and  use  coding.  Another  method 
for  simple  messages  is  to  use  an  algorithm  which  constructs  a message 
out  of  correctly  received  segments.  The  results  of  this  work  is  in 
terms  of  message  reception  probabilities.  Message  reception  probabilities 
are  ultimately  the  most  desired  result  and  the  algorithms  to  calculate 
these  results  are  being  constructed. 
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SECTION  VII 


I 

SUMMARY 


ttU  report  has  demonstrated  methods  to  calculate  message  reception 
probabilities  and  other  useful  infection  starting  from  the  statistical 
description  of  the  propagation  environment.  A discussion  on  the  state 
of  knowledge  of  the  ambient  environments  is  also  given.  The  presented 
methods  provide  sufficient  accuracy  to  examine  the  various  quest 
affecting  the  satellite  —cations  systems  mentioned  in  the  int.o 

duction. 
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ALTITUDE  (Km) 

Figure  7.  Ionospheric  Parameters 
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TIME  ABOVE  ORDINATE 


Figure  12.  FSK  Link  Performance 


MAXIMUM  WINDOW 


Figure  13.  PSK  Link  Performance 


Appendix  A 

Generation  of  Complex  Correlated  Random  Sequences 


Let  f(t)  be  a zero  mean  random  Gaussian  variable  on  the  domain 
(-T,T).  f(t)  can  be  expanded  in  a Fourier  series. 


£(t)-|S  cn  - f1 

n=-°° 


1 f “iWnt 

=n  = f-  J f(t)  A n dt 


G(t-t  ) = f*(t)f(t  ) 


.J(tTt').  = fCt)f(t') 


LimCTt^  Cm)  = 0,  m^n,  m ^ 0 


LimCTC^  C m)  - 4 I da  G(a)  cos  (au>m) 


= T K 


Lim(TCn  Cm)  =0,ra^n,  m^O 


Lim(TCm  C_ 


t-  1 

-J  ' « i 


'n  da  F(a)  cos  (aw  ) 
U m 


- + i v 


If  T is  chosen  such  that  F(t)  = G(t)  ss  0 for  t •>  T,  the  limits  in 

equations  A5,  A6,  A7,  and  A8  can  be  assumed.  Now  if  C ■ C + iC  . , and 

m mr  mi 

using  equations  A5,  A6,  A7  and  A8,  then  with  m t 0. 


C C = C . C 4 = K /2 
mr  mr  mi  mi  m' 


CL,,.  C . = 0 
mr  ^mi 


C C = -C  . C = L /2 
mr  -mr  mi  -m1  m' 


C C , = C Cml  = M /2 
mr  -mi  -mr  mi  _ 


Since  C , C C , and  C . are  random  Gaussian  variables,  the 
mr  mi  -mr  -mi 

joint  probability  distribution  function  of  these  variables  can  be 

written  in  teems  of  , and  M..  .....  ■ ...  . ... 

m m m 


^^mr ’^mi’^-mr’^*^i^  ~ ^ exP 


C C . C C . 
mr  mi  -mr  “mi 


This  expression  can  be  simplified  by  diagonalizing  the  4x4  matrix. 
To  diagonalize  the  above  matrix,  a transformation  matrix  must  be  found 
with  the  following  properties . 


1 


where  I is  the  identity  matrix,  Q is  the  4x4  matrix  in  equation  13, 


and  R is  a diagonal  matrix  with  diagonal  elements  A^,  A2,  A^,  and  A^.  T 


diagonal  elements  are  easily  solved  for. 


h ' X2  ‘ Km  + ^ + ^ 


A matrix  T,  with  the  desired  properties  is 


M //“ 
m 


m 


-L  hT 
m 


-*Jr 


-V^ 


where 


^ 2 

= + . The  probability  density  function  can 


be  written  as 


, » ft  I 

P(Cm  , cm  , c , c mi) 

v m r*  m i’  -mr  ~mi 


= C exp 


L p<c'mr2  + C'-I2)  + 

L^Sn  + *V  + “m2 


2(C ' 2 + c'  2) 

-mr  -mi 


**  - K2  + «m 


' C -1  > ~ 
+ M?  . 


J 


mm 


where 


c 

mr 

1 

0 

-K^ 

-v^ 

c„, 

mi 

= 1 

0 

1 

-M  /✓“ 
m 

C-mr 

/2 

l tr 

m 

m tr 

m 

1 

0 

C-mi 

\'r 

0 

1 

I 

Cmr 

Cmi 

V 

^ -mr 

C’-mi 

The  primed  variables  are  also  random  Gaussian  variables  with  variances 
as  shown  in  equation  A19.  Generating  a simple  of  f(t)  on  the  interval 

-T  < t < T is  now  straight  forward.  For  each  m t 0,  sample  a Gaussian 

» * » » 

distribution  four  times  to  get  values  of  C mr>  C mi»  C i an<*  ^ -mi’ 

transform  using  equation  20  to  unprimed  variables,  and  form  the  random 

Fourier  coefficient  by  C = Cmr  + iC  . . For  m = 0,  another  procedure 

m mr  mi 

must  be  used.  From  equations  A6  and  A8 


^or  c0i  ~ Mo/2 

Thus 


+ V'2 

‘ <K0  - V'2 
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(A20) 


(A21) 

(A22) 

(A23) 

(A24) 

(A25) 


The  joint  probability  function  can  now  be  written  as 
p(Cor>  Col>  = C exP  j[lColcc 


Jor  i 


(K0-L0) 

-Mo 


-M. 


(Ko+Lo) 


c 

oi 

C 

or 

as 

with  m 

/(k2-l2-m2) 

O 0 o 


be  diagonalized. 
Let 


1/2 


cor  * [(VLo)/21"  (Cor  + C'ol) 


o"0 


Col  - («0-Lo)/2I1/2  - c'oj) 


Then 


p(Cor»  Goi)  = c exP 


| - l [ 2(Ko  - LP  - "o  «o  - Lo>1/2)  c'or 


2 2 9 
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(A27) 
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Appendix  B 


Noise  Derived  Power  Spectral  Density  of  a Square  Law  Device 

Let  Y(t)  = a(s(t)  + n(t))2 

where  s(t)  and  n(t)  are  the  signal  and  noise  input  respectively  to  the 
square  law  device. 

Y(t)  = a[s2(t)  + 2n(t)s(t)  + n2(t)] 

If  n(t)  is  a zero  mean, stationary,  gaussian  process  and  independent  of 

s(t)  and  if  s(t)  is  stationary,  then 
«* 

Y = a[s2  + n2] 

Also 

Y = a[s  + 6 s4  n + n ] 

The  autocorrelation  function  of  Y(t)  is 

RyOO  = a2[(s(t)  + n(t))2(s(t  + T)  + n(t  + T))2] 

= a2[s2(t)s2(t  + T)  + n2(t)n^(t  + x) 

+ 4 s(t)s(t  + x)  n(t)n(t  + X)  + 2 s2(t)  n2(t)] 
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Ry(T)  = az[Rg2(x)  + 2 s2"  n2  + 4 ^(1)^(1)  + R^x)  ] 


The  last  three  terms  in  equation  B7  are  noise  derived.  The  power 
spectral  density  of  Y(t)  is  the  Fourier  transform  of  Py(x) . The 
spectral  density  associated  with  the  middle  terms  is 


Ssx„  (£)  ' 4a2/  V£> 


. -i2nfx 


+ 2a2  s2  n^  6(f) 


Is „ <£' 


) ss(f-f')  df' 


+ 2a2  s2  n2  6(f) 


Since  n(t)  is  Gaussian, 


-Rn2  (x)  = 2a2  R^  (x)  + a2  s' 


Thus  the  spectral  density  of  the  last  term  is 


S (f)  = 2a2  [ S (f ')  Sn(f-f’)  df' 
nxn  J n 


2 T2 

+ a s2 


The  total  power  spectral  density  from  noise  related  terms  for  the  square 


law  device  is 


SsXn<£)  + Snxn(£>  ' 43  S"  <£)  ‘ S*  (£) 


+ 2a2S_  (£)  * S.(f)  + [s2  + s2  n2]a26(f) 


where  denotes  a convolution  operation.  It  is  important  to  note  that 
the  process  with  a spectral  density  of  (f ) * Sn (f ) is  not  gaussian  as 
is  the  Sn(f)  *S  (f)  process. 


Appendix  C 

SOLUTION  TECHNIQUES  FOR  THE  FIRST  AND  SECOND  ORDER  LOOP  ERROR  EQUATIONS 
The  equation  for  the  second  order  loop  equation  to  be  solved  is 


de(Q  = d6Jt!  _ KNJfm  gin  (e(t))  _ JL  n(t) 


dt 


dt 


- as(  ~ sin  (e(u))  + 

! r\ 


n(u) 

\jf7i 


du 


(Cl) 


where  n(t)  is  a zero  mean  Gaussian  process  with  a power  spectral  density 
of  N/2  centered  at  f = 0.  Let 


F(t) 


(C2) 


Then 

F(t+At)  = F(t)  - aK^~  sin  (e(t))At  - Ti  (t+At)At  (C3) 

° VPo 

j,  ft+At 

where  I (t+At)  t2  = \ n(u)du. 

Equation  1 can  be  similarly  integrated  for  a result  for  e(t+At) 
where  it  is  assumed  that  At  is  sufficiently  small  so  that  P(t)  and 
sin(e(t))  do  not  change  significantly  over  At. 


e(t+At)  = e(t)  + 


d6CO  _ K>/Ipl  sin  ( e(t) ) + F (t) 


dt 


At 


. aK^JUtT  ln  (c(t))At2  - I (t+At)At3/2  - -pb  Ix(t+At )At!<  (C4) 
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2 P 


54 


55 


For  e(t+At) 


e(t+At)  = e(t)  + 


where 


/ 2 


a 


2 


I2(t+At) 


a„  = 


+ F (t)  + 


a 


4 


c 

2 


The  overall  solution  method  is  now  a straightforward  procedure 
First  the  initial  conditions  are: 


F(0)  = 0 
e(0)  = 0 

Next  pick  an  initial  At.  Then  proceed  through  the  following. 

1.  Calculate  and 

2.  Solve  equation  BIO 

3.  Calculate  anc*  a4* 


4.  Solve  equation  B13 


5.  Check  size  of  time  step. 

If  | e(t+At)  -e(t)  | > 20°,  then  At  = 3At/4  and 

If  | e(t+At)  -e(t)  | < 5°,  then  At  = 4At/3,  increment  total  time, 
and  go  to  step  1. 

The  method  of  solution  outlined,  while  somewhat  inelegant,  does  pro- 
vide a reasonable  simulation  of  the  nonlinear  phase  lock  loop  behavior  in 
the  presence  of  noise  while  remaining  simple  and  efficient  on  a digital 
computer. 

The  method  of  solution  of  the  first  order  loop  is  the  same  as  for 
the  second  order  loop  with  a = 0. 
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